pis_robotics_control_systemfandomcom-20200215-history
Simulating Physical System with Simulink
This is a preliminary guide to simulate a dynamical system using Simulink. The set-up of the code is similar. However, we will use different solver to find solutions to our system. The one of benefits of Simulink is the intuitive interface , which allows the user to connect multiple system together without further derivation. We will look at the first example, an RC low-pass filter. The same set of questions also apply. There are several things to keep in mind for writing the code for simulation. # How long does the simulation is going to take ? # What is the smallest time step of the simulation ? Are they fixed ? # Do you have the state-space representation of your system ? # What is the initial condition of your system ? # Is your system a closed-loop system ? any inputs ? # How will you interpret your result ? animation, or graph ? Simulink Environment and Workspace First, we will use MATLAB to store information. Then, we will use Simulink to simulate the dynamics as well as export the result back to the workspace. We will generate 1 MATLAB file as well 1 Simulink Model. The format of the code that I use is the following. # Setting up parameters and initial condition in MATLAB # Defining the dynamics of the system in Simulink # Call the simulation in MATLAB # Interpret the solution with graph and animation in MATLAB These are criteria for the simulation of the RC filter. * The initial time is zero. The duration of simulation is 10 seconds. * We are interested in observing the voltage across the capacitor as well as the current through the circuit. * Initially, the switch is off for very long time. * The switch is turned on right at the beginning. * The power supply generate a sinusoidal input with the natural frequency of 1 Hz and the amplitude of 5 V. * The resistance is 10Ω. The capacitance of the circuit is 1mF. * The power supply, resistor, and capacitor are connected in series. * The dynamics of the RC circuit can be described with the following state-space representation. \dot{x_1} = \frac{u -x_1}{\tau} where u =A\sin(\omega t) \tau =RC x_1 is the voltage across the capacitor in V u is the voltage from power supply in V R is resistance of the resistor in \Omega C is capacitance of the capacitor in F It might not be intuitive, but in order to find solutions with Simulink, the differential equation has to be transformed into an integral equation. The integral equation also requires one to specify the initial condition of the system as well. \int_{t_0}^t \! \dot{x_1} \, \mathrm{d}t. = \int_{x_1(t_0)}^{x_1(t)} \! \, \mathrm{d}x_1. = x_1(t)-x_1(t_0) = \int_{t_0}^t \! \frac{u-x_1}{\tau} \, \mathrm{d}t. By substituting initial conditions, the state variable x_1 can be expressed in an integral form. x_1(t) = \int_{t_0}^t \! \frac{u-x_1}{\tau} \, \mathrm{d}t.+v_0 (1) where v_0 is the initial voltage across the capacitor Setting Up Parameters and Initial Condition in MATLAB In case of testing the simulation multiple times with different value of parameters, it is more convenient to declare variables in the beginning of the script. These variables will store in MATLAB's workspace. One can access, change, and clear these variable at any time in the code as long as they are declared before. For the sake of consistency, we will define parameters into 2 different categories. # Physical Parameters # Initial Condition Different from the previous approach, the time element will be adjusted in Simulink. Physical Parameters Remember to use the right units. The following is the sample code in MATLAB for declaring physical parameters of the system. R = 10; C = 1*10^-3; A = 5; w = 2*pi*1; tau = R*C; Initial Conditions The initial condition is 0 since the capacitor was not charged prior to the initial time. The following is the sample code in MATLAB for assigning the initial conditions of the system. We will denote the initial condition by "IC". IC = 0; Defining the dynamics of the system in Simulink Internal Dynamics Create a new SImulink model. A Simulink workspace will appear. We are going to implement the equation (1) with variety of blocks. Before adding anymore, "Model Configuration Parameters" has to be checked and adjust. "Model Configuration Parameters" can be found under "Simulation" tab. We will leave as it is for now. Notice that the start time and end time are defaulted as 0 and 10, respectively. Model Configuration Parameters The first block that we will use is the "Integrator" block. This can be found in "Commonly Used Blocks". Place the block in the workspace and double-click the block. Change the initial condition from 0 to IC. Notice: to be explicit, you can also make the initial condition external. Function Block Parameters: Integrator At this point, we will generate a mathematic expression from blocks, which is equivalent to the integrand in equation (1). Notice: the expression inside of the integral is the rate of change of voltage across the capacitor. expression = \dot{x_1} = \frac{u-x_1}{\tau} Summation block or "Sum" block can also be found in "Commonly Used Blocks". Add a sum block and double-click it. You will see two operators defaulted as "++" in the pop-up window. Change them to "+-" and click OK. Also add a gain block to the workspace and double-click it. Adjust the gain to 1/tau and click OK. Connect the output of the sum block to the input of the gain block. Then connect the output of the gain block with the input of the integrator block. At this point, you will have 3 blocks connected in series. Notice: you can name the signal to prevent confusion. Cascaded Block To add a feedback x_1 , we simply connect the output of the integrator block back to the negative input of the sum block. Closed-loop Internal Dynamics Since we are interested in the voltage across the capacitor, we can attach an "output" block to the end of the integrator block. Closed-loop Internal Dynamics with Output Masking Parameters This step is not necessary for our case. However, masking parameters is useful for duplicating and re-adjusting sub-system. Before masking the parameters, select all components except the output x_1. Press Ctrl+G in order to make the selected components a sub-system. Double-click the sub-system and rename both input and output to "u" and "x_1", respectively. Notice: you can now navigate through different level of the system via the left-hand-side "Model Browser". Subsystem with Model Browser Subsystem with Input & Output We will also rename the sub-system as RC filter. Navigate to the highest level and right-click on the RC filter. Select Mask>Create Mask. Select "Parameters" tab and add a parameter. Set the prompt and variable to "Time Constant" and "tau", respectively. Click Ok. Mask Editor: RC Filter Now you can double click RC Filter. A function block parameter will appear. You can put "tau" as the parameter as well as "R*C". Creating an input Sometimes, the input of a system is time-dependent. In our case, the input of the system is a sinusoidal wave that depend on time. u = A\sin(\omega t) You can either type "clock" in the search bar or find it under "sources". "Clock" will generate time signal for the entire simulation. Add a clock and another gain block to the workspace. Double-click on the gain block and change the gain to "w". Connect the clock to the input of the gain block. Find "Trigonometric Function" block under "Math Operations". Add the trig function block to the workspace. Also make sure that function is set to "sine". Connect the output of the gain block to the input of the trig function block. Then, add another gain block at the end of the trig function block. Adjust the gain of that block to "A". At this point, you should have 4 blocks connected in series and a sub-system. Unconnected Input and Sub-system In order to simplify the model, we will transform 4 cascaded blocks into a sub-system called "Sinusoidal Input". We will also mask "A" and "w" to "Amplitude" and "Angular Frequency". Now, we can connected two sub-system together. Two Sub-systems Obtaining results using scope and export Now that we have a model to simulate, it is important to know what are we looking for. As stated in the beginning that we are interested in observing both voltage across the capacitor and the current through the circuit, it is necessary to determine these quantities. It is obvious that the voltage across the capacitor is x_1. However, the current through the circuit have different behavior. Since all components are connected in series, the current through all components are the same. Therefore, we can say that the current through the capacitor is equal to the current through the entire circuit. Hence, we can express the current as the following. i_c = C\frac{dv_c}{dt} where i_c is the current through the capacitor. One might think to add a "derivative" block with x_1. However, this can take additional processing power. Therefore, an alternative way to find derivative of x_1 is to obtain the expression before the integrator block. We will go back to RC Filter and modify the block. Notice: you can no longer double-click to access the block. Right-click on RC Filter and select "Open in New Tab". Add another gain block to the workspace. Double-click the gain block and change the gain to "C". Connect the input of the gain block directly to the line in front of the integrator block. Also add another output block at the end of the gain block and rename it to "i_c". Notice: now, sub-system RC Filter has 2 outputs at the highest level of the model. RC Filter with 2 outputs RC Filter with input and 2 outputs We can combine two results to make the output more compact by using "Mux". Add a mux to the workspace. Connect both input of the mux to the output of the filter. At this point, we can observe the result using scope or store them in MATLAB workspace. If you decide to use a scope, simply add a scope block to the workspace and connect it with the output of the mux. You must run the script that contains the physical parameters before running the simulation. An alternative way to obtain the data is to use "To Workspace" block. This block can be found in "Sinks". Add "To Workspace" block to the workspace and double-click the block to change the name of the variable. Once the simulation is run, all information of x_1 and i_c will be stored in the variable that specified in "To Workspace" block. Save the model as "RCFilter". The complete model of RC Circuit Call the simulation in MATLAB You do not have to press "Run" in the Simulink interface all the time. You can run the Simulink model using a MATLAB command called "sim". Given the name of the model as an arguement, "sim" will return a Simulink signal, which is a structure that consists of various data structure. In order to get a vector/ matrix of data, you can use "." followed by the name of the field. For example, result.Data will return nx2 matrix, where n is the amount of the sample data. The following is the sample code in MATLAB for simulating and obtaining data from a Simulink Model. sim('RCFilter'); t = result.time; x_1 = result.Data(:,1); i_c = result.Data(:,2); Interpret the Solutions with Graph and Animation For this example, we will learn how to create an animation of moving graph rather than static graph. The idea behind the animation is to display one frame after the other in sequence. In MATLAB, we will set the first frame, which includes a regular static plot . After that, we will use for-loop to iterate through each data point and change the setting of the original plot. A function "pause" is required in order to make the animation run according to the time vector. First, we plot 1 set of data. After that, we increase the amount of data by a fixed number of step. Then, we re-plot again. The following is a sample code for full simulation of RC Filter and animation of the moving graph. R = 10; C = 1*10^-3; A = 5; w = 2*pi*1; tau = R*C; IC = 0; sim('RCFilter'); t = result.time; x_1 = result.Data(:,1); i_c = result.Data(:,2); step = 1; xprevious = x_1(1); xcurrent = x_1(1+step); iprevious = i_c(1); icurrent = i_c(1+step); figure; hold on; xdata = t(1:1); x_1data = x_1(1:1); i_cdata = i_c(1:1); pX = plot(xdata,x_1data,'b'); pI = plot(xdata,i_cdata,'r'); axis(10 -5 5); for i = 1+step:step:length(t), xdata = t(1:i); x_1data = x_1(1:i); i_cdata = i_c(1:i); hold on; set(pX,'XData',xdata,'YData',x_1data); set(pI,'XData',xdata,'YData',i_cdata); pause(t(i)-t(i-step)); end